//
// Created by sunmy on 2020/10/17.
//

#ifndef HW1_CHOLESKY_DECOMPOSITION_H
#define HW1_CHOLESKY_DECOMPOSITION_H

#include "Eigen/Core"
#include <iostream>

template <int n = 10, typename T = long double>
void cholesky_decomposition_impl()
{
    using namespace Eigen;
    using namespace std;
    typedef Matrix<T, n, n> MatrixND;
    typedef Matrix<T, n, 1> VectorND;

    MatrixND H = MatrixND();

    for (auto i = 0; i < n; ++i)
    {
        for (auto j = 0; j < n; ++j)
        {
            H(i, j) = (T)1.0 / (i + j + 1);
        }
    }

    VectorND b = H * VectorND::Ones(n, 1);

    //    cout << "Now the matrix H is:" << endl << H << endl;
    //    cout << "Now the vector b is:" << endl << b << endl;

    MatrixND L = MatrixND::Zero(n, n);

    for (auto j = 0; j < n; ++j)
    {
        T sum = H(j, j);
        for (auto k = 0; k < j; ++k)
            sum -= L(j, k) * L(j, k);
        L(j, j) = sqrt(sum);
        for (auto i = j + 1; i < n; ++i)
        {
            T s = H(i, j);
            for (auto k = 0; k < j; ++k)
                s -= L(i, k) * L(j, k);
            L(i, j) = s / L(j, j);
        }
    }

    VectorND y = VectorND();

    for (auto i = 0; i < n; ++i)
    {
        y(i) = b(i);
        for (auto j = 0; j < i; ++j)
            y(i) -= y(j) * L(i, j);
        y(i) /= L(i, i);
    }

    VectorND x = VectorND();
    for (auto i = n - 1; i >= 0; --i)
    {
        x(i) = y(i);
        for (auto j = i + 1; j < n; ++j)
            x(i) -= x(j) * L(j, i);
        x(i) /= L(i, i);
    }

    cout.precision(12);
    cout << "The solution for n = " << n << " is:" << endl
         << x << endl;
}

void cholesky_decomposition()
{
    cholesky_decomposition_impl<10>();
    cholesky_decomposition_impl<11>();
    cholesky_decomposition_impl<12>();
    cholesky_decomposition_impl<13>();
    cholesky_decomposition_impl<14>();
    cholesky_decomposition_impl<15>();
}

#endif //HW1_CHOLESKY_DECOMPOSITION_H
